This lab uses the prioritize R package for systematic conservation planning which explores some of the common tasks involved with developing prioritizations.

if (!require("librarian")){
  install.packages("librarian")
  library(librarian)
}
librarian::shelf(
  assertthat, BiocManager, dplyr, gridExtra, here, mapview, 
  prioritizr, prioritizrdata, 
  raster, remotes, rgeos, rgdal, scales, sf, sp, stringr,
  units)
if (!require("lpsymphony")){
  BiocManager::install("lpsymphony")
  library(lpsymphony)
}

Data Setup

dir_data <- here("data/prioritizr")
pu_shp <- file.path(dir_data, "pu.shp")
pu_url <- "https://github.com/prioritizr/massey-workshop/raw/main/data.zip"
pu_zip <- file.path(dir_data, basename(pu_url))
vegetation_tif <- file.path(dir_data, "vegetation.tif")

dir.create(dir_data, showWarnings = F, recursive = T)
if(!file.exists(pu_shp)){
  download.file(pu_url, pu_zip)
  unzip(pu_zip, exdir = dir_data)
  dir_unzip <- file.path(dir_data, "data")
  files_unzip <- list.files(dir_unzip, full.names = T)
  file.rename(
    files_unzip,
    files_unzip %>% str_replace("prioritizr/data", "prioritizr"))
  unlink(c(pu_zip, dir_unzip), recursive = T)
}

Data

Data Import

The data contains vector-based planning unit data (pu.sh) and raster-based data describing the spatial distribution of 32 vegetation classes (vegetation.tif) in southern Tasmania, Australia.

# import planning unit data
pu_data <- as(read_sf(pu_shp), "Spatial")

# format columns in planning unit data
pu_data$locked_in <- as.logical(pu_data$locked_in)
pu_data$locked_out <- as.logical(pu_data$locked_out)

# import vegetation data
veg_data <- stack(vegetation_tif)

Planning unit data

The planning unit data contains spatial data describing the geometry for each planning unit and attribute data including: - id: unique identifiers for each planning unit - cost: acquisition cost values for each planning unit (millions of Australian dollars) - status: status information for each planning unit (only relevant with Marxan) - locked_in: logical values (ie. TRUE/FALSE) indicating if planning units are covered by protected areas or not - locked_out: logical values (ie. TRUE/FALSE) indicating if planning units cannot be managed as a protected area because they are too degraded.

# print a short summary of the data
print(pu_data)
## class       : SpatialPolygonsDataFrame 
## features    : 516 
## extent      : 348703.2, 611932.4, 5167775, 5344516  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=55 +south +datum=WGS84 +units=m +no_defs 
## variables   : 5
## names       :   id,             cost, status, locked_in, locked_out 
## min values  :  557, 3.59717531470679,      0,         0,          0 
## max values  : 1130,  47.238336402701,      2,         1,          1
# plot the planning unit data
plot(pu_data)

# plot an interactive map of the planning unit data
mapview(pu_data)
# print the structure of object
str(pu_data, max.level = 2)
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  516 obs. of  5 variables:
##   ..@ polygons   :List of 516
##   ..@ plotOrder  : int [1:516] 69 104 1 122 157 190 4 221 17 140 ...
##   ..@ bbox       : num [1:2, 1:2] 348703 5167775 611932 5344516
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
# print the class of the object
class(pu_data)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
# print the slots of the object
slotNames(pu_data)
## [1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"
# print the coordinate reference system
print(pu_data@proj4string)
## CRS arguments:
##  +proj=utm +zone=55 +south +datum=WGS84 +units=m +no_defs
# print the first six rows of the data
head(pu_data@data)
##    id     cost status locked_in locked_out
## 1 557 29.74225      0     FALSE      FALSE
## 2 558 29.87703      0     FALSE      FALSE
## 3 574 28.60687      0     FALSE      FALSE
## 4 575 30.83416      0     FALSE      FALSE
## 5 576 38.75511      0     FALSE      FALSE
## 6 577 38.11618      2      TRUE      FALSE
# print the first six values in the cost column of the attribute data
head(pu_data$cost)
## [1] 29.74225 29.87703 28.60687 30.83416 38.75511 38.11618
# print the highest cost value
max_pu_cost <- max(pu_data$cost)
max_pu_cost
## [1] 47.23834
# print the smallest cost value
min(pu_data$cost)
## [1] 3.597175
# print the average cost value
mean(pu_data$cost)
## [1] 26.87393
# plot a map of the planning unit cost data
spplot(pu_data, "cost")

# plot an interactive map of the planning unit cost data
mapview(pu_data, zcol = "cost")

Question 1: How many planning units are in the planning unit data?

num_pu <- nrow(pu_data@data)

Answer 1: The planning unit data contains 516 planning units.

Question 2: What is the highest cost value? Answer 2: The highest cost planning units has a values of $47.2383364 million Austrailian dollars.

Question 3: Is there a spatial pattern in the planning unit cost values (hint: use plot to make a map)?

Answer 3: Yes. The lower cost planning units are located in the north east section of the study area. The most expensive planning units are located in the north central and north west sections of the study area.

Vegetation data

The vegetation data is in raster formate and describes the spatial distribution of 32 vegetation classes in the study area. The raster data contains multiple layers (also called “bands”) and each layer corresponds to a spatial grid with exactly the same area adn has exactly the same dimensionality (ie number of rows, columns and cells). In this dataset there are 32 different regular spatial grids layered on top of each other with each layer corresponding to a different vegetation class. Each layer contains a grin with 164 rows, 326 columns, and 53464 cells. Within each layer, each cell corresponds to a 0.967 by 1.02 km square. The values associated with each grid cell indicate the (1) presence or (0) absence of a given vegetation class in the cell.

# print a short summary of the data
print(veg_data)
## class      : RasterStack 
## dimensions : 164, 326, 53464, 32  (nrow, ncol, ncell, nlayers)
## resolution : 967, 1020  (x, y)
## extent     : 298636.7, 613878.7, 5167756, 5335036  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=55 +south +datum=WGS84 +units=m +no_defs 
## names      : vegetation.1, vegetation.2, vegetation.3, vegetation.4, vegetation.5, vegetation.6, vegetation.7, vegetation.8, vegetation.9, vegetation.10, vegetation.11, vegetation.12, vegetation.13, vegetation.14, vegetation.15, ... 
## min values :            0,            0,            0,            0,            0,            0,            0,            0,            0,             0,             0,             0,             0,             0,             0, ... 
## max values :            1,            1,            1,            1,            1,            1,            1,            1,            1,             1,             1,             1,             1,             1,             1, ...
# plot a map of the 20th vegetation class
plot(veg_data[[20]])

# plot an interactive map of the 20th vegetation class
mapview(veg_data[[20]])
# print number of rows in the data
nrow(veg_data)
## [1] 164
# print number of columns in the data
ncol(veg_data)
## [1] 326
# print number of cells in the data
ncell(veg_data)
## [1] 53464
# print number of layers in the data
nlayers(veg_data)
## [1] 32
# print resolution on the x-axis
xres(veg_data)
## [1] 967
# print resolution on the y-axis
yres(veg_data)
## [1] 1020
# print spatial extent of the grid, ie coordinates for corners
extent(veg_data)
## class      : Extent 
## xmin       : 298636.7 
## xmax       : 613878.7 
## ymin       : 5167756 
## ymax       : 5335036
# print the coordinate reference system
print(veg_data@crs)
## CRS arguments:
##  +proj=utm +zone=55 +south +datum=WGS84 +units=m +no_defs
# print a summary of the first layer in the stack
print(veg_data[[1]])
## class      : RasterLayer 
## band       : 1  (of  32  bands)
## dimensions : 164, 326, 53464  (nrow, ncol, ncell)
## resolution : 967, 1020  (x, y)
## extent     : 298636.7, 613878.7, 5167756, 5335036  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=55 +south +datum=WGS84 +units=m +no_defs 
## source     : vegetation.tif 
## names      : vegetation.1 
## values     : 0, 1  (min, max)
# print the value in the 800th cell in the first layer of the stack
print(veg_data[[1]][800])
##   
## 0
# print the value of the cell located in the 30th row and the 60th column of the first layer
print(veg_data[[1]][30, 60])
##   
## 0
# calculate the sum of all the cell values in the first layer
cellStats(veg_data[[1]], 'sum')
## [1] 17
# calculate the maximum value of all the cell values in the first layer
cellStats(veg_data[[1]], 'max')
## [1] 1
# calculate the minimum value of all the cell values in the first layer
cellStats(veg_data[[1]], 'min')
## [1] 0
# calculate the mean value of all the cell values in the first layer
cellStats(veg_data[[1]], 'mean')
## [1] 0.00035239

Question 1: What part of the study area is the 13th vegetation class found (hint: make a map)? For instance, is it in the south-east part of the study area?

# plot a map of the 13th vegetation class
plot(veg_data[[13]])

Answer 1: The 13th vegetation class is mostly found in the north east section of the study area.

Question 2: What proportion of cells contain the 12th vegetation class?

plot(veg_data[[12]])

mapview(veg_data[[12]])
cells_with_12th_class <- (cellStats(veg_data[[12]], 'sum') / ncell(veg_data[[12]])) * 100
cells_with_12th_class
## [1] 1.531872

Answer 2: Approximately 1.53% of cells contain the 12th vegetation class.

Question 3: Which vegetation class is the most abundant (ie present in the greatest number of cells)?

num_layers <- nlayers(veg_data)
veg_layers <- seq(from = 1, to = num_layers, by = 1)
veg_class_abun <- data.frame()
for (i in seq_along(veg_layers)) {
  abun <- cellStats(veg_data[[i]], 'sum')
  abun_df <- data.frame(abun)
  veg_class_abun <- rbind(veg_class_abun, abun_df)
}
veg_class_most_abun <- which.max(veg_class_abun$abun)
veg_class_most_abun
## [1] 12
num_cell_12th <- cellStats(veg_data[[veg_class_most_abun]], "sum")

Answer 3: The 12th vegetation class is the most abundant and is present in the greatest number of cells compared to the other vegetation layers. There are 819 cells with the 12th vegetation class present.

# cellStats(veg_data, "sum")

Gap analysis

Introduction

This gap analysis involves calculating how well each biodiversity feature is represented (covered) by protected areas. The project compares current representation of each feature by protected areas (defined as 5% of spatial distribution covered by protected areas) to a target threshold (defined as 20% of their spatial distribution covered by protected areas). The target threshold denotes the minimum amount (ie minimum proportion of spatial distribution) that we need of each feature to be represented in the protected area system.

Feature abundance

Calculate how much of each vegetation feature occurs inside each planning unit (ie the abundance of the feature).

# create prioritizr problem with only the data
p0 <- problem(pu_data, veg_data, cost_column = "cost")

# print empty problem
# we can see that only the cost and feature data are defined
print(p0)
## Conservation Problem
##   planning units: SpatialPolygonsDataFrame (516 units)
##   cost:           min: 3.59718, max: 47.23834
##   features:       vegetation.1, vegetation.2, vegetation.3, ... (32 features)
##   objective:      none
##   targets:        none
##   decisions:      default
##   constraints:    <none>
##   penalties:      <none>
##   portfolio:      default
##   solver:         default
# calculate amount of each feature in each planning unit
abundance_data <- feature_abundances(p0)

# print abundance data
print(abundance_data)
## # A tibble: 32 × 3
##    feature       absolute_abundance relative_abundance
##    <chr>                      <dbl>              <dbl>
##  1 vegetation.1                16.0                  1
##  2 vegetation.2                14.3                  1
##  3 vegetation.3                10.4                  1
##  4 vegetation.4                17.8                  1
##  5 vegetation.5                13.0                  1
##  6 vegetation.6                14.3                  1
##  7 vegetation.7                20.0                  1
##  8 vegetation.8                14.0                  1
##  9 vegetation.9                18.0                  1
## 10 vegetation.10               20.0                  1
## # … with 22 more rows
# note that only the first ten rows are printed
# this is because the abundance_data object is a tibble (ie tbl_df) object and not a standard data.frame object
print(class(abundance_data))
## [1] "tbl_df"     "tbl"        "data.frame"
# we can print all of the rows in abundance_data like this
print(abundance_data, n = Inf)
## # A tibble: 32 × 3
##    feature       absolute_abundance relative_abundance
##    <chr>                      <dbl>              <dbl>
##  1 vegetation.1                16.0                  1
##  2 vegetation.2                14.3                  1
##  3 vegetation.3                10.4                  1
##  4 vegetation.4                17.8                  1
##  5 vegetation.5                13.0                  1
##  6 vegetation.6                14.3                  1
##  7 vegetation.7                20.0                  1
##  8 vegetation.8                14.0                  1
##  9 vegetation.9                18.0                  1
## 10 vegetation.10               20.0                  1
## 11 vegetation.11               23.6                  1
## 12 vegetation.12              748.                   1
## 13 vegetation.13              126.                   1
## 14 vegetation.14               10.5                  1
## 15 vegetation.15               17.5                  1
## 16 vegetation.16               15.0                  1
## 17 vegetation.17              213.                   1
## 18 vegetation.18               14.3                  1
## 19 vegetation.19               17.1                  1
## 20 vegetation.20               21.4                  1
## 21 vegetation.21               18.6                  1
## 22 vegetation.22              297.                   1
## 23 vegetation.23               20.3                  1
## 24 vegetation.24              165.                   1
## 25 vegetation.25              716.                   1
## 26 vegetation.26               24.0                  1
## 27 vegetation.27               18.8                  1
## 28 vegetation.28               17.5                  1
## 29 vegetation.29               24.7                  1
## 30 vegetation.30               59.0                  1
## 31 vegetation.31               60.0                  1
## 32 vegetation.32               32                    1

The absolute_abundance column contains the total amount of each feature in all the planning units. The relative_abundance column contains the total amount of each feature in the planning units expressed as a proportion of the total amount in the underlying raster data. Since all the raster cells containing vegetation overlap with the planning units, all of the values in the relative_abundance column are equal to one (meaning 100%).

# add a new column with the feature abundances expressed in area units (ie km^2)
abundance_data$absolute_abundance_km <- (abundance_data$absolute_abundance * prod(res(veg_data))) %>% 
  set_units(m^2) %>% 
  set_units(km^2)

# print abundance data
print(abundance_data)
## # A tibble: 32 × 4
##    feature       absolute_abundance relative_abundance absolute_abundance_km
##    <chr>                      <dbl>              <dbl>                [km^2]
##  1 vegetation.1                16.0                  1                  15.8
##  2 vegetation.2                14.3                  1                  14.1
##  3 vegetation.3                10.4                  1                  10.2
##  4 vegetation.4                17.8                  1                  17.6
##  5 vegetation.5                13.0                  1                  12.8
##  6 vegetation.6                14.3                  1                  14.1
##  7 vegetation.7                20.0                  1                  19.7
##  8 vegetation.8                14.0                  1                  13.9
##  9 vegetation.9                18.0                  1                  17.8
## 10 vegetation.10               20.0                  1                  19.7
## # … with 22 more rows

Explore the abundance data

# calculate the average abundance of the features
mean(abundance_data$absolute_abundance_km)
## 86.82948 [km^2]
# plot a histogram of the feature abundances
hist(abundance_data$absolute_abundance_km, main = 'Feature Abundances')

# find the name of the feature with the largest abundance
abundance_data$feature[which.max(abundance_data$absolute_abundance_km)]
## [1] "vegetation.12"

Question 1:What is the median abundance of the features (hint median)?

median_abundance <- median(abundance_data$absolute_abundance_km)
median_abundance
## 19.1165 [km^2]

Answer 1: The median abundance of the features is 19.12km^2.

Question 2: What is the name of the feature with smallest abundance?

min_abun_feature <- abundance_data$feature[which.min(abundance_data$absolute_abundance_km)]
min_abun_feature
## [1] "vegetation.3"

Answer 2: The feature with the smallest abundance is vegetation.3.

Question 3: How many features have a total abundance greater than 100 km^2 (hint: use sum(abundance_data$absolute_abundance_km2 > set_units(threshold, km^2) with the correct threshold value)

num_features_greater100km <- sum(abundance_data$absolute_abundance_km > set_units(100, km^2))
num_features_greater100km
## [1] 6

Answer 3: There are 6 features that have a total abundance greater than 100 km/^2.

Feature Representation

After calculating the total amount of each feature in the planning units (ie the features abundance), calculate the amount of each feature in the planning units that are covered by protected ares (ie feature representation by protected areas) using the eval_feature_representation_summary() function. This function requires (i) a conservation problem object with the planning unit and biodiversity data and also (ii) an object representing a solution to the problem (i.e an object in the same format as the planning unit data with values indicating if the planning units are selected or not).

# create column in planning unit data wtih binary values (zeros and ones) indicating if a planning unit is covered by protected areas or not
pu_data$pa_status <- as.numeric(pu_data$locked_in)

# calculate feature representation by protected areas
repr_data <- eval_feature_representation_summary(p0, pu_data[, "pa_status"])

# print feature representation data
print(repr_data)
## # A tibble: 32 × 5
##    summary feature       total_amount absolute_held relative_held
##    <chr>   <chr>                <dbl>         <dbl>         <dbl>
##  1 overall vegetation.1          16.0         0            0     
##  2 overall vegetation.2          14.3         0            0     
##  3 overall vegetation.3          10.4         0            0     
##  4 overall vegetation.4          17.8         0            0     
##  5 overall vegetation.5          13.0         0            0     
##  6 overall vegetation.6          14.3         0            0     
##  7 overall vegetation.7          20.0         0            0     
##  8 overall vegetation.8          14.0         0            0     
##  9 overall vegetation.9          18.0         0.846        0.0470
## 10 overall vegetation.10         20.0         0            0     
## # … with 22 more rows

The absolute_held column shows the total amount of each feature held in the solution (ie the planning units covered by protected areas). The relative_held column shows the proportion of each feature held in the solution (ie the proportion of each feature’s spatial distribution held in protected areas). Since the absolute_held values correspond to the number of grid cells in the veg_data object with overlap with protected areas, convert then to area units (ie km^2) so they can be reported.

# add new column with the areas represented in km^2
repr_data$absolute_held_km <- (repr_data$absolute_held * prod(res(veg_data))) %>% 
  set_units(m^2) %>% 
  set_units(km^2)

# print representation data
print(repr_data)
## # A tibble: 32 × 6
##    summary feature       total_amount absolute_held relative_held absolute_held_km
##    <chr>   <chr>                <dbl>         <dbl>         <dbl>           [km^2]
##  1 overall vegetation.1          16.0         0            0                 0    
##  2 overall vegetation.2          14.3         0            0                 0    
##  3 overall vegetation.3          10.4         0            0                 0    
##  4 overall vegetation.4          17.8         0            0                 0    
##  5 overall vegetation.5          13.0         0            0                 0    
##  6 overall vegetation.6          14.3         0            0                 0    
##  7 overall vegetation.7          20.0         0            0                 0    
##  8 overall vegetation.8          14.0         0            0                 0    
##  9 overall vegetation.9          18.0         0.846        0.0470            0.834
## 10 overall vegetation.10         20.0         0            0                 0    
## # … with 22 more rows

Investigate how well the species are represented.

Question 1: What is the average proportion of the features held in protected area (hint: use mean(table$relative_held) with the correct table name)?

avg_prop_protected <- (mean(repr_data$relative_held)) * 100
avg_prop_protected
## [1] 24.15487

Answer 1: An average of 24.155% of features are protected.

Question 2: If we set a target of 10% coverage by protected areas, how many features fail to meet this target (hint: use sum(table$relative_held >= target_vale) with the correct table name)?

num_features_lessthan_10pct <- sum(repr_data$relative_held <= 0.1)
num_features_lessthan_10pct
## [1] 17

Answer 2: Based on current protected status, 17 features fail to meet a target of 10% coverage by protected areas.

Question 3: If we set a target of 20% coverage by protected areas, how many features fail to meet this target?

num_features_lessthan_20pct <- sum(repr_data$relative_held <= 0.2)
num_features_lessthan_20pct
## [1] 18

Answer 3: Based on current protected status, 18 features fail to meet a target of 20% coverage by protected areas.

Question 4: Is there a relationship between the total abundance of a feature and how well it is represented by protected areas?

hint: plot(abundance_data$absolute_abundance ~ repr_data$relative_held)

plot(abundance_data$absolute_abundance ~ repr_data$relative_held)

Answer 4: There does not appear to be a relationship between the total abundance of a feature and how well it is represented by protected areas.

Spatial Prioritizations

Introduction

Develop prioritizations to identify priority areas for protected area establishment. Note, prioritizr is a decision support tool similar to Marxan and Zonation. It is designed to help you make decisions - it can’t make decisions for you.

Starting out simple

Create a prioritization using the minimum set formulation of the reserve selection problem. This formulation means that we want a solution that will meet the targets for our biodiversity features for minimum cost. Here, we will set 5% targets for each vegetation class and use the data in the cost column to specify acquisition costs.

# print planning unit data
print(pu_data)
## class       : SpatialPolygonsDataFrame 
## features    : 516 
## extent      : 348703.2, 611932.4, 5167775, 5344516  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=55 +south +datum=WGS84 +units=m +no_defs 
## variables   : 6
## names       :   id,             cost, status, locked_in, locked_out, pa_status 
## min values  :  557, 3.59717531470679,      0,         0,          0,         0 
## max values  : 1130,  47.238336402701,      2,         1,          1,         1
# make prioritization problem
p1_rds <- file.path(dir_data, "p1.rds")
if(!file.exists(p1_rds)){
  p1 <- problem(pu_data, veg_data, cost_column = "cost") %>% 
    add_min_set_objective() %>% 
    add_relative_targets(0.05) %>% # representation targets
    add_binary_decisions() %>% 
    add_lpsymphony_solver()
  saveRDS(p1, p1_rds)
}

p1 <- readRDS(p1_rds)

# print problem
print(p1)
## Conservation Problem
##   planning units: SpatialPolygonsDataFrame (516 units)
##   cost:           min: 3.59718, max: 47.23834
##   features:       vegetation.1, vegetation.2, vegetation.3, ... (32 features)
##   objective:      Minimum set objective 
##   targets:        Relative targets [targets (min: 0.05, max: 0.05)]
##   decisions:      Binary decision 
##   constraints:    <none>
##   penalties:      <none>
##   portfolio:      default
##   solver:         Lpsymphony [first_feasible (0), gap (0.1), time_limit (2147483647), verbose (1)]
# solve problem
s1 <- solve(p1)

# print solution, the solution_1 column contains the solution values
# indicating if a planning unit is (1) selected or (0) not
print(s1)
## class       : SpatialPolygonsDataFrame 
## features    : 516 
## extent      : 348703.2, 611932.4, 5167775, 5344516  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=55 +south +datum=WGS84 +units=m +no_defs 
## variables   : 7
## names       :   id,             cost, status, locked_in, locked_out, pa_status, solution_1 
## min values  :  557, 3.59717531470679,      0,         0,          0,         0,          0 
## max values  : 1130,  47.238336402701,      2,         1,          1,         1,          1
# calculate number of planning units selected in the prioritization
eval_n_summary(p1, s1[, "solution_1"])
## # A tibble: 1 × 2
##   summary  cost
##   <chr>   <dbl>
## 1 overall    15
# calculate total cost of the prioritization
eval_cost_summary(p1, s1[, "solution_1"])
## # A tibble: 1 × 2
##   summary  cost
##   <chr>   <dbl>
## 1 overall  385.
# plot solution
# selected = green, not selected = grey
spplot(s1, "solution_1", col.regions = c("grey80", "darkgreen"), main = "s1, colorkey = FALSE")

Examine the solution.

Question 1: How many planning units were selected in the prioritization? What proportion of planning units were selected in the prioritization?

num_pu_selected_1 <- as.numeric(eval_n_summary(p1, s1[, "solution_1"]))
num_pu_selected_1 <- num_pu_selected_1[2]

prop_pu_selected_1 <- (num_pu_selected_1 / num_pu) * 100
prop_pu_selected_1
## [1] 2.906977

Answer 1: A total of 15 planning units were selected in the prioritization. This represents 2.907% of all planning units.

Question 2: Is there a pattern in the spatial distribution of the priority areas?

Answer 2: There does not appear to be a pattern in the spatial distribution of the priority areas. The planning units selected in the prioritization are scattered across the study area.

Question 3: Can you verify that all of the targets were met in the prioritization? hint: `eval_feature_representation_summary(p1, s1[, “solution_1”])

s1_summary <- eval_feature_representation_summary(p1, s1[, "solution_1"])
min_relative_held_s1 <- (min(s1_summary$relative_held)) * 100
min_relative_held_s1
## [1] 5.331146

Answer 3: The minimum of a feature held in the solution is 5.331% so all targets were met in the prioritization.

Adding complexity

Add additional constraints to the problem to make it more useful. First, lock in planning units that are already covered by protected areas. If some vegetation communities are already secured inside existing protected areas, then we might not need to add as many new protected areas to the existing protected area system to meet targets.

# plot locked_in data
# TRUE = blue, FALSE = grey
spplot(pu_data, "locked_in", col.regions = c("grey80", "darkblue"), main = "locked_in", colorkey = FALSE)

# make prioritization problem
p2_rds <- file.path(dir_data, "p2.rds")
if (!file.exists(p2_rds)){
  p2 <- problem(pu_data, veg_data, cost_column = "cost") %>% 
    add_min_set_objective() %>% 
    add_relative_targets(0.05) %>% 
    add_locked_in_constraints("locked_in") %>% 
    add_binary_decisions() %>% 
    add_lpsymphony_solver()
  saveRDS(p2, p2_rds)
}
p2 <- readRDS(p2_rds)

# print problem
print (p2)
## Conservation Problem
##   planning units: SpatialPolygonsDataFrame (516 units)
##   cost:           min: 3.59718, max: 47.23834
##   features:       vegetation.1, vegetation.2, vegetation.3, ... (32 features)
##   objective:      Minimum set objective 
##   targets:        Relative targets [targets (min: 0.05, max: 0.05)]
##   decisions:      Binary decision 
##   constraints:    <Locked in planning units [198 locked units]>
##   penalties:      <none>
##   portfolio:      default
##   solver:         Lpsymphony [first_feasible (0), gap (0.1), time_limit (2147483647), verbose (1)]
# solve problem
s2 <- solve(p2)

# print solution
print(s2)
## class       : SpatialPolygonsDataFrame 
## features    : 516 
## extent      : 348703.2, 611932.4, 5167775, 5344516  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=55 +south +datum=WGS84 +units=m +no_defs 
## variables   : 7
## names       :   id,             cost, status, locked_in, locked_out, pa_status, solution_1 
## min values  :  557, 3.59717531470679,      0,         0,          0,         0,          0 
## max values  : 1130,  47.238336402701,      2,         1,          1,         1,          1
# plot solution
# selected = green, not selected = grey
spplot(s2, "solution_1", col.regions = c("grey80", "darkgreen"), main = "s2, colorkey = FALSE")

Now set the target to 10%.

# make prioritization problem
p3_rds <- file.path(dir_data, "p3.rds")
if(!file.exists(p3_rds)){
  p3 <- problem(pu_data, veg_data, cost_column = "cost") %>% 
    add_min_set_objective() %>% 
    add_relative_targets(0.1) %>% 
    add_locked_in_constraints("locked_in") %>% 
    add_binary_decisions() %>% 
    add_lpsymphony_solver()
  saveRDS(p3, p3_rds)
}
p3 <- readRDS(p3_rds)

# print problem
print(p3)
## Conservation Problem
##   planning units: SpatialPolygonsDataFrame (516 units)
##   cost:           min: 3.59718, max: 47.23834
##   features:       vegetation.1, vegetation.2, vegetation.3, ... (32 features)
##   objective:      Minimum set objective 
##   targets:        Relative targets [targets (min: 0.1, max: 0.1)]
##   decisions:      Binary decision 
##   constraints:    <Locked in planning units [198 locked units]>
##   penalties:      <none>
##   portfolio:      default
##   solver:         Lpsymphony [first_feasible (0), gap (0.1), time_limit (2147483647), verbose (1)]
# solve problem
s3 <- solve(p3)

# print solution
print(s3)
## class       : SpatialPolygonsDataFrame 
## features    : 516 
## extent      : 348703.2, 611932.4, 5167775, 5344516  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=55 +south +datum=WGS84 +units=m +no_defs 
## variables   : 7
## names       :   id,             cost, status, locked_in, locked_out, pa_status, solution_1 
## min values  :  557, 3.59717531470679,      0,         0,          0,         0,          0 
## max values  : 1130,  47.238336402701,      2,         1,          1,         1,          1
# plot solution
# selected = green, not selected = grey
spplot(s3, "solution_1", col.regions = c("grey80", "darkgreen"), main = s3, colorkey = FALSE)

Include locked_out planning areas.

# plot locked_out data
# TRUE = red, FALSE = grey
spplot(pu_data, "locked_out", col.regions = c("grey80", "darkred"), main = "locked_out", colorkey = FALSE)

# make prioritization problem
p4_rds <- file.path(dir_data, "p4.rds")
if(!file.exists(p4_rds)){
  p4 <- problem(pu_data, veg_data, cost_column = "cost") %>% 
    add_min_set_objective() %>% 
    add_relative_targets(0.1) %>% 
    add_locked_in_constraints("locked_in") %>% 
    add_locked_out_constraints("locked_out") %>% 
    add_binary_decisions() %>% 
    add_lpsymphony_solver()
  saveRDS(p4, p4_rds)
}
p4 <- readRDS(p4_rds)

# print problem
print(p4)
## Conservation Problem
##   planning units: SpatialPolygonsDataFrame (516 units)
##   cost:           min: 3.59718, max: 47.23834
##   features:       vegetation.1, vegetation.2, vegetation.3, ... (32 features)
##   objective:      Minimum set objective 
##   targets:        Relative targets [targets (min: 0.1, max: 0.1)]
##   decisions:      Binary decision 
##   constraints:    <Locked in planning units [198 locked units]
##                    Locked out planning units [6 locked units]>
##   penalties:      <none>
##   portfolio:      default
##   solver:         Lpsymphony [first_feasible (0), gap (0.1), time_limit (2147483647), verbose (1)]
# solve problem
s4 <- solve(p4)

# print solution
print(s4)
## class       : SpatialPolygonsDataFrame 
## features    : 516 
## extent      : 348703.2, 611932.4, 5167775, 5344516  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=55 +south +datum=WGS84 +units=m +no_defs 
## variables   : 7
## names       :   id,             cost, status, locked_in, locked_out, pa_status, solution_1 
## min values  :  557, 3.59717531470679,      0,         0,          0,         0,          0 
## max values  : 1130,  47.238336402701,      2,         1,          1,         1,          1
# plot solution
# selected = green, not selected = grey
spplot(s4, "solution_1", col.regions = c("grey80", "darkgreen"), main = s4, colorkey = FALSE)

Compare the solutions

Question 1: What is the cost of the planning units selected in s2, s3, and s4?

cost_s2 <- eval_cost_summary(p2, s2[, "solution_1"])
cost_s2 <- cost_s2[2]

cost_s3 <- eval_cost_summary(p3, s3[, "solution_1"])
cost_s3 <- cost_s3[2]

cost_s4 <- eval_cost_summary(p4, s4[, "solution_1"])
cost_s4 <- cost_s4[2]

Answer 1: The cost of the selected planning units is 6600.09 million Australian dollars for s2, 6669.91 million Australian dollars for s3, and 6711.58 million Australian dollars for s4.

Question 2: How many planning units are in s2, s3, and s4?

num_pu_s2 <- eval_n_summary(p2, s2[, "solution_1"])
num_pu_s2 <- num_pu_s2[2]

num_pu_s3 <- eval_n_summary(p3, s3[, "solution_1"])
num_pu_s3 <- num_pu_s3[2]

num_pu_s4 <- eval_n_summary(p4, s4[, "solution_1"])
num_pu_s4 <- num_pu_s4[2]

Answer 2: There are 205 planning units in s2, 211 planning units in s3, and 212 planning units in s4.

Question 3: Do the solutions with more planning units have a greater cost? Why (or why not)?

Answer 3: The solutions with more planning units have a greater cost. This is because, by excluding the locked out areas from the solution, more planning units were required to achive conservation targets. Also, the locked at areas my have lower values since they are degraded.

Question 4: Why does the first solution (s1) cost less than the second solution with protected areas locked into the solution (s2)?

Answer 4: The protected areas could have higher cost values. The agorithm is setup to minimuze costs, but when the low cost degraded planning units are locked out, conservation targets must be met with higher cost areas.

Question 5: Why does the third solution (s3) cost less than the fourth solution with the highly degraded areas locked out (s4)?

Answer 5: The highly degraded areas likely have lower acquisition cost values. By not being able to include these planning units in the s4 solution, the fourth solution costs more than the third solution.

Penalizing fragmentation

To promote connectivity within the planned protected area system, add penalties to the conservation planning problem to penalize fragmentation. Specify a trade-off between the primary objective (solution cost) and fragmentation (total exposed boundary length) using a penalty value. Generally, use penalty values between 0.00001 and 0.01. Since the planning unit data is in a spatial format (ie. vector or raster data), prioritizr can automatically calculate the boundary data.

# make prioritization problem
p5_rds <- file.path(dir_data, "p5.rds")
if(!file.exists(p5_rds)){
  p5 <- problem(pu_data, veg_data, cost_column = "cost") %>% 
    add_min_set_objective() %>% 
    add_boundary_penalties(penalty = 0.001) %>% 
    add_relative_targets(0.1) %>% 
    add_locked_in_constraints("locked_in") %>% 
    add_locked_out_constraints("locked_out") %>%
    add_binary_decisions() %>% 
    add_lpsymphony_solver()
  saveRDS(p5, p5_rds)
}
p5 <- readRDS(p5_rds)

# print problem
print(p5)
## Conservation Problem
##   planning units: SpatialPolygonsDataFrame (516 units)
##   cost:           min: 3.59718, max: 47.23834
##   features:       vegetation.1, vegetation.2, vegetation.3, ... (32 features)
##   objective:      Minimum set objective 
##   targets:        Relative targets [targets (min: 0.1, max: 0.1)]
##   decisions:      Binary decision 
##   constraints:    <Locked out planning units [6 locked units]
##                    Locked in planning units [198 locked units]>
##   penalties:      <Boundary penalties [edge factor (min: 0.5, max: 0.5), penalty (0.001), zones]>
##   portfolio:      default
##   solver:         Lpsymphony [first_feasible (0), gap (0.1), time_limit (2147483647), verbose (1)]
# solve problem
# note: this will take longer to run than the previous runs
s5 <- solve(p5)
# print solution
print(s5)
## class       : SpatialPolygonsDataFrame 
## features    : 516 
## extent      : 348703.2, 611932.4, 5167775, 5344516  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=55 +south +datum=WGS84 +units=m +no_defs 
## variables   : 7
## names       :   id,             cost, status, locked_in, locked_out, pa_status, solution_1 
## min values  :  557, 3.59717531470679,      0,         0,          0,         0,          0 
## max values  : 1130,  47.238336402701,      2,         1,          1,         1,          1
# plot solution
# selected = green, not selected = grey
spplot(s5, "solution_1", col.regions = c("grey80", "darkgreen"), main = s5, colorkey = FALSE)

Compare the solutions to the problems with (s5) and without (s4) the boundary length penalties.

Question 1: What is the cost of the fourth (s4) and fifth (s5) solutions? Why does the fifth solution (s5) cost more than the fourth (s4) solution?

cost_s4 <- eval_cost_summary(p4, s4[, "solution_1"])
cost_s4 <- cost_s4[2]

cost_s5 <- eval_cost_summary(p5, s5[, "solution_1"])
cost_s5 <- cost_s5[2]
num_pu_s4 <- eval_n_summary(p4, s4[, "solution_1"])
num_pu_s4 <- num_pu_s4[2]

num_pu_s5 <- eval_n_summary(p5, s5[, "solution_1"])
num_pu_s5 <- num_pu_s5[2]

Answer 1: Solution 4 costs 6711.58 million Australian dollors and solution 5 costs 6747.59 million Australian dollars. Solution 5 costs more due to the boundary length penalty. The planning units that minimize the boundary length cost more than other planning units that would achieve the target vegetation goals. Solution 5 also includes more planning units than solution 4.

Question 2: Try setting the penalty value to 0.000000001 (i.e. 1e-9) instead of 0.001. What is the cost of the solution now? Is it different from the fourth solution (s4)? Hint: try plotting the solutions to visualize them. Is this a useful penalty value? Why (or why not)?

# make prioritization problem
p6_rds <- file.path(dir_data, "p6.rds")
if(!file.exists(p6_rds)){
  p6 <- problem(pu_data, veg_data, cost_column = "cost") %>% 
    add_min_set_objective() %>% 
    add_boundary_penalties(penalty = 0.000000001) %>% 
    add_relative_targets(0.1) %>% 
    add_locked_in_constraints("locked_in") %>% 
    add_locked_out_constraints("locked_out") %>%
    add_binary_decisions() %>% 
    add_lpsymphony_solver()
  saveRDS(p6, p6_rds)
}
p6 <- readRDS(p6_rds)

# print problem
print(p6)
## Conservation Problem
##   planning units: SpatialPolygonsDataFrame (516 units)
##   cost:           min: 3.59718, max: 47.23834
##   features:       vegetation.1, vegetation.2, vegetation.3, ... (32 features)
##   objective:      Minimum set objective 
##   targets:        Relative targets [targets (min: 0.1, max: 0.1)]
##   decisions:      Binary decision 
##   constraints:    <Locked in planning units [198 locked units]
##                    Locked out planning units [6 locked units]>
##   penalties:      <Boundary penalties [edge factor (min: 0.5, max: 0.5), penalty (1e-09), zones]>
##   portfolio:      default
##   solver:         Lpsymphony [first_feasible (0), gap (0.1), time_limit (2147483647), verbose (1)]
s6 <- solve(p6)
# plot solution
# selected = green, not selected = grey
spplot(s6, "solution_1", col.regions = c("grey80", "darkgreen"), main = s6, colorkey = FALSE)

cost_s6 <- eval_cost_summary(p6, s6[, "solution_1"])
cost_s6 <- cost_s6[2]

num_pu_s6 <- eval_n_summary(p6, s6[, "solution_1"])
num_pu_s6 <- num_pu_s6[2]

Answer 2: With a lower boundary penalty, the cost of the solution is 6711.76 million Australian dollars and the solution includes 211 planning units. This solution has a slightly higher cost than the fourth solution. The fourth solution and solution with the lower boundary penalty both have the same number of planning units. Visually, they appear very similar with only a few differences in selected planning units. This penalty value is not useful because there are still several isolated planning units that are not connected to other units identified for conservation.

Question 3: Try setting the penalty value to 0.5. What is the cost of the solution now? Is it different from the fourth solution (s4)? Hint: try plotting the solutions to visualize them? Is this a useful penalty value? Why (or why not)?

# make prioritization problem
p7_rds <- file.path(dir_data, "p7.rds")
if(!file.exists(p7_rds)){
  p7 <- problem(pu_data, veg_data, cost_column = "cost") %>% 
    add_min_set_objective() %>% 
    add_boundary_penalties(penalty = 0.5) %>% 
    add_relative_targets(0.1) %>% 
    add_locked_in_constraints("locked_in") %>% 
    add_locked_out_constraints("locked_out") %>%
    add_binary_decisions() %>% 
    add_lpsymphony_solver()
  saveRDS(p7, p7_rds)
}
p7 <- readRDS(p7_rds)

# print problem
print(p7)
## Conservation Problem
##   planning units: SpatialPolygonsDataFrame (516 units)
##   cost:           min: 3.59718, max: 47.23834
##   features:       vegetation.1, vegetation.2, vegetation.3, ... (32 features)
##   objective:      Minimum set objective 
##   targets:        Relative targets [targets (min: 0.1, max: 0.1)]
##   decisions:      Binary decision 
##   constraints:    <Locked in planning units [198 locked units]
##                    Locked out planning units [6 locked units]>
##   penalties:      <Boundary penalties [edge factor (min: 0.5, max: 0.5), penalty (0.5), zones]>
##   portfolio:      default
##   solver:         Lpsymphony [first_feasible (0), gap (0.1), time_limit (2147483647), verbose (1)]
s7 <- solve(p7)
# plot solution
# selected = green, not selected = grey
spplot(s7, "solution_1", col.regions = c("grey80", "darkgreen"), main = s7, colorkey = FALSE)

cost_s7 <- eval_cost_summary(p7, s7[, "solution_1"])
cost_s7 <- cost_s7[2]

num_pu_s7 <- eval_n_summary(p7, s7[, "solution_1"])
num_pu_s7 <- num_pu_s7[2]

Answer 3: The solution with a boundary length penalty of 0.5 has a cost of 9816.64 million Australian dollars and includes 314 planning units. This solution is very different from the fouth solution because now all selected planning units are connected (there are no isolated selected units). Since a single conservation area this large and with such a high cost is likely not feasible, the 0.5 penalty value is not useful from a feasibility perspective. It is only useful if reducing fragmentation is a top priority.